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A LOW COMPLEXITY ALGORITHM FOR NON-MONOTONICALLY 

EVOLVING FRONTS * 

ALEXANDRA TCHENG, JEAN-CHRISTOPHE NAVE t 

Abstract. A new algorithm is proposed to describe the propagation of fronts advected in the 
normal direction with prescribed speed function F. The assumptions on F are that it does not 
depend on the front itself, but can depend on space and time. Moreover, it can vanish and change 
sign. To solve this problem the Level-Set Method [Osher, Sethian; 1988] is widely used, and the 
Generalized Fast Marching Method [Carlini et al.; 2008] has recently been introduced. The novelty 
of our method is that its overall computational complexity is predicted to be comparable to that 
of the Fast Marching Method [Sethian; 1996], [Vladimirsky; 2006] in most instances. This latter 
algorithm is 0(N n log N n ) if the computational domain comprises N n points. Our strategy is to use 
it in regions where the speed is bounded away from zero - and switch to a different formalism when 
F ~ 0. To this end, a collection of so-called sideways partial differential equations is introduced. 
Their solutions locally describe the evolving front and depend on both space and time. The well- 
posedness of those equations, as well as their geometric properties are adressed. We then propose 
a convergent and stable discretization of those PDEs. Those alternative representations are used to 
augment the standard Fast Marching Method. The resulting algorithm is presented together with 
a thorough discussion of its features. The accuracy of the scheme is tested when F depends on 
both space and time. Each example yields an 0(1/N) global truncation error. We conclude with a 
discussion of the advantages and limitations of our method. 

Key words, front propagation, Hamilton-Jacobi equations, fast marching method, level-set 
method, optimal control, viscosity solutions. 
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1. Introduction. The design of robust numerical schemes describing front prop¬ 
agation has been a subject of active research for several decades. The need for such 
schemes is felt across many areas of applied sciences: geometric optics [36], optimal 
control [25, 53], lithography [2, 3, 4], shape recognition [31, 29], dendritic growth 
[26, 45], gas and fluid dynamics [32, 33, 52], combustion [58], etc. Depending on 
the problem at hand, various issues may arise. Consider the following two interface 
propagation phenomena: A fire propagating through a forest, and a large evolving 
population of bacteria in a Petri dish. In either case, space can be divided into distinct 
regions: burnt vs. unburnt, and populated vs. unpopulated. The boundaries between 
those regions form fronts that evolve in time. Those examples differ from one another 
in that a fire front can only propagate monotonically , whereas bacteria may advance 
or recede, depending on the stimuli present in their environment. This distinction 
led to different approaches when modelling those evolutions. Monotone propaga¬ 
tion can be recast into a ‘static’ problem, as opposed to non-monotone evolution, 
which is instrinsically time-dependent. As a result, efficient single-pass algorithms 
for monotone propagation have been developed. In contrast, accurate algorithms for 
non-monotonically evolving fronts require a larger number of computations. In this 
paper, we propose a model that reconciles the advantages of previous methods - We 
accurately describe non-monotone front evolution with an algorithm that performs a 
low number of operations. 

One of the early means of accurately propagating fronts was to use the Level-Set 
Method (LSM) [35]. This implicit approach embeds the front as the zero-level-set of 
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an auxiliary function <fi. In the above example, <f> could be negative in regions occupied 
by bacteria, and positive in other regions. Each contour of this level-set function is 
then evolved under the given speed function F, which guarantees that the front itself 
moves properly. The robustness and simplicity of the first order discretization of 
this problem made it popular. Additionally, this approach can handle a very wide 
class of speed functions, including those that change sign. However, describing the 
evolution of an (n — l)-dimensional front in 1" requires solving for a function of 
n + 1 variables, since <f> depends on space as well as time. Moreover, in order for 
the solution to remain accurate, it is often desirable to enforce the signed distance 
property |V(/>| « 1 in a neighbourhood of the front. There exists a vast literature 
on lowering the computational complexity of the LSM, cf. [5, 39, 44, 37], and on 
maintaining the accuracy of the solution, cf. [11, 12, 39, 51, 40]. Nevertheless, those 
features are incorporated at the expense of the simplicity and the efficiency of the 
original LSM. 

The Fast Marching Method (FMM) [41, 54] constitutes the second significant 
advance in the field. This approach requires the speed function to be bounded away 
from zero, and to be only space-dependent. Under those conditions, the FMM builds 
the ‘first arrival time’ function i/i such that to every point x in space is associated the 
value t = Y’(x) at which the front reaches x, cf. [41, 42, 47, 44], In the context of 
fire propagation, ijj records the time at which the parcel of land burnt. The use of a 
Dijkstra-like data structure [19] renders this scheme very efficient. A variant of this 
algorithm known as the Fast Sweeping Method runs in 0(N n ) complexity [57] when 
the computational domain comprises N n points. Recently, Falcone et al. [9] proposed 
a Generalized FMM (GFMM) that is able to handle vanishing speeds. This algorithm 
is supported by theoretical results on its convergence in the class of viscosity solutions. 
The examples presented are found to accurately propagate the fronts subject to a wide 
range of speed functions. However, when F depends on time, the GFMM no longer 
makes use of a Dijkstra-like data structure. Its overall complexity is expected to revert 
to that of the LSM in such instances. 

In the light of this previous work, it is desirable to design an algorithm able to 
handle speed functions that change sign, while retaining the efficiency of the FMM. 
This is the main purpose of this article. Note that if F changes sign, a point x in 
space may be reached by the front several times. This implies that the arrival time 
can no longer be described as a function depending solely on space. However, it 
is still possible to locally describe it as the graph of a function. Consider the set 
Af := {(x, t) : x belongs to the front at time t}. The set A4 consists of the surface 
traced out by the fronts as they evolve through space and time. If Ad embeds as a 
C fc -manifold of dimension n in K ra x (0,T), then by definition, each point (x, t) £ Af 
belongs to a neighbourhood that is locally the image of a C fc -function of n variables. 
The fact that under mild assumptions A4 is a compact subset of x (0, T) guarantees 
that we only need a finite number of neighbourhoods to cover Ad, or equivalently, a 
finite number of functions to parametrize Ad. The images of those functions - which 
possibly depend on time as well as space - provide local representations of the set Ad . 
Our approach makes use of those other representations whenever the purely spatial 
one is not available - e.g., when n = 2 and Ad cannot be locally described by the 
standard first arrival time function {t = ip(x, y)}, we may describe it as {x = <)} 

or {y = ^(a:,t)}. To this end, we introduce sideways PDEs solved by those C k - 
functions. We illustrate in detail how they relate to previous work, argue that they 
are well-posed, and show that their solution does provide a local description of AL 
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Moreover, we provide a scheme to discretize them, prove that it converges to the 
correct viscosity solution, and show that it is stable. 

In practice, the proposed algorithm amounts to augmenting the FMM to be able 
to describe M near those points (x, t) where F(x, t) = 0. The fact that different 
representations are used to build different parts of M implies that those pieces need 
to be woven together along their overlapping parts, to form a single codimension 
one subset of R” x (0,T). This is done by storing the (n + l)-dimensional normal 
associated to each point and by using interpolation. To illustrate the overall method, 
examples are presented where an 0(1/N) global truncation error is achieved. Those 
tests all feature speed functions that vanish, and possibly depend on time. 

Since the algorithm always approximates a function of n-variables, the dimen¬ 
sionality of the problem is never raised, unlike what happens in the LSM. As a result, 
the computational complexity is expected to be comparable to that of the FMM. 

Outline of the article. This paper is organized as follows. We state the problem 
we are addressing in §2. We also present the LSM and the FMM, before providing 
a simple example to motivate our method. The case where F is bounded away from 
zero and depends on time is addressed in §3. The sideways PDEs we use in regions 
where F ps 0 are introduced in §4. A discussion of their properties is provided along 
with a convergent and stable scheme to discretize them. We explain how the different 
formalisms can be woven into a single method in §5. The pseudo-codes are given and 
discussed in § 6 . We predict the complexity and accuracy of the overall method in 
§7 and § 8 . Four examples are then covered in details in § 8 . Those assess the global 
behaviour and the accuracy of the scheme. The advantages and weaknesses of our 
approach are discussed in §9, where an additional example is covered to address the 
limitations of the method. We conclude in §10. 

2. Preliminaries. 

2.1. Problem statement. Let a subset Co C 1" be closed with no boundary. 
Assume it is an orientable manifold of codimension one, with a well defined unique 
outer normal fio(x). Suppose Co is advected in time, and denote the resulting subset 
of R” at time t by Ct . We want to describe Ct for 0 < t < T in the case where each 
point x G C t is advected under the velocity 


v = u(x, t) = F(x, f)n(x, t) (2.1) 

i.e., with the prescribed speed function F = F(x, f), in the direction of the outward 
normal to Ct, n = n(x,£). 

2.2. Assumptions. In addition to the assumptions already stated, in the rest 
of this paper we assume that the following hold. The initial set Co is known exactly, 
and is assumed to be C 2 in the sense that if it is given as the image of a map, 
e.g., 7 : 5 n_1 —> R ra , then 7 G C 2 (5' ra_1 ). The speed F = F(x, t) is known exactly 
for all (x, t). Unless otherwise specified, it is allowed to vanish and change sign. It 
does not depend on the curve itself, or any of its derivatives. For simplicity, we also 
make the following strong assumption: the map F : R n x (0, T) —> R is analytic. In 
particular, this implies that the subset defined as F := {(x, t) : F(x, t) = 0 } is closed 
and has codimension one in R n x [0,T]. We let K be the Lipschitz constant of F. 
Together, those assumptions guarantee that for any given t G (0, T), there exists a 
well defined normal n = n(x, t ) almost everywhere along Cj. 
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2.3. Previous Work. For completeness we briefly go over two of the meth¬ 
ods mentioned in the introduction. Considering that the set R” \ C t consists of two 
connected components, we define At to be the bounded one. 

2.3.1. The Level-Set Method. This approach was introduced by Osher & 
Sethian in [35]. Their idea is to embed the curve Ct as the zero-level-set of a function 4> : 
R" x [0, T] —s> R, i.e., C t = {x : </>(x, t) = 0}. In this setting, the outward normal n(x, t) 
is The Level-Set Equation is derived from linear advection (f> t + F(x, t) ■ V0 = 0 

to yield the following Initial Value Problem (IVP): 


( fa+F |V0| = 0 on R ra x (0, T) 

\ </>(x, 0) = </>o(x) on R” x {0} 

where </>o(x) is such that {x : c/>o(x) = 0} = Cq. This method enjoys many desirable 
properties that have been studied in a variety of contexts [20, 21, 22, 23, 37, 44]. One 
of the most prominent is that topological changes are accurately handled, and do 
not require special treatment. In [35], the authors propose various discretizations of 
this evolution on a spatial domain that comprises N n points. The resulting method 
has complexity 0(N n ) at each time step, due to the fact that all the contours of the 
level-set function are advected. To lower this high computational cost, it is possible to 
work only within a neighbourhood of the zero-level-set: This yields the Narrow Band 
LSM [5]. To be able to render the curve Ct accurately, it is desirable to preserve the 
signed distance property | V0| ~ 1. To this end, the reinitialization method has been 
studied extensively [39, 40, 51, 11]. Early versions of this method tend to displace the 
zero-level-set, yielding inaccuracies in the final Ct- Moreover, they usually involve a 
large number of computations. 

2.3.2. The Fast Marching Method. The Fast Marching Method was inde¬ 
pendently proposed by Sethian [41] & Tsitsiklis [54]. Strongly rooted in control theory, 
it requires that F = F(x) > <5 > 0 on R". Under those conditions, the FMM solves 
the following Eikonal equation, whose unknown is the time ^ : R" i->Rat which each 
point is reached by the curve 


f l v ^l = j? °n Aq \ Cq 

\ ?/>(x) = 0 on Cq 


The FMM makes use of a Narrow Band to advance the front in a manner that enforces 
the characteristic structure of the PDE into the solution. See [41, 42, 47, 44] and [25] 
for details. Recent improvements of this method include on the one hand the work of 
Zhao [57], who further lowered the complexity of the algorithm to develop the Fast 
Sweeping Method. On the other hand, Vladimirsky relaxed the restrictions on the 
speed by allowing it to be time-dependent. We discuss this latter method in §3. 

2.4. Motivation. We first present a simple example to motivate the need for 
an augmented FMM. Consider the initial curve Cq = {x : x 2 + y 2 = Tq} C R 2 and 
the time-dependent speed F(t) = 1 — ct, where c and ro are positive constants. Let 
</>o(x) be the signed distance function (/>o(x) = ^/x 2 + y 2 — ro =: r(x) — r 0 . The 
exact solution to the IVP (2.2) is then cj)(x,y,t) = r(5c) — (r 0 — (ct 2 /2 — t)). The 
evolution of the curve can be formally split into two parts: (1) For t £ [0, ^], the 
circle expands until it reaches the maximal radius R = ro + J-C- ( 2 ) Fort £ (i,T], 
where T = — - (l — yT + 2cro), the circle contracts until it collapses to the point 
(0, 0) at time T. 
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Fig. 2.1. Chart decomposition of Ai when c — 2 and ro = 0.25. The circle collapses to (0,0) 
at time T fa 1.21. (0) The manifold, sliced by the plane t = 0.83 to yield the magenta curve Co. 83 - 
Three other typical curves Ct are featured with dashed lines. Co appears as a thick plain line. (1) 
Ai along with Wi,— and Wi,+ appearing in green. (2) Ai along with W 2 ,— and W 2 ,+ appearing in 
blue. (3) Ai along with W 3 ,— and W 3 ,+ appearing in red. 


Consider the following atlas si to describe the resulting C 0 -manifold Ai featured 
on Figure 2.1. Let W:=lx [0,T]. Then si — U?_ 1 {('0i ) ±, Wj,±)} where the real¬ 
valued functions '>pi t ± are defined as: 


and 


V’i- 

: U — 

> [~R,o] 

fpi,+ 

u — 

■>[0 ,R] 

1p2- 

: U — 

> [-i?,0] 

*p2,+ 

u — 


1p3- 

: R 2 - 

-+ [o, k) 

V’3,+ 

R 2 - 



(2.4) 


= ±\J(r 0 - ct 2 /2 + tf - y 2 (2.5) 

ip 2 ,±(x,t) = ±\J(ro - ct 2 /2 + t) 2 - x 2 (2.6) 

ip3,±{x,y) = \ ( x ± \A ~ 2 c(r(x) - r 0 )) (2.7) 


We also define the sets W,.± as the real part of the image of the functions Those 
sets are featured on Figure 2.1. The functions ipa t ± can be verihed to be the unique 
classical solutions to: 

= F& 3,-(g)) 0n ^3 

= 0 on Co 


|VV>3,-(x)| 

^3,-(x) 


( 2 . 8 ) 
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r iw 3 ,+(x)| = onU ^+ f 2 9 ) 

1 tfe,+(x) = \ on C\/ c 

where U 3 - = {x : tq < r(x) < i?} and Uz i+ = {x : 0 < r(x) < i?}. Together, the 
graphs of ip%- and ^> 3 ,+ describe all of Ad but the circle of radius R reached at time 
t = -. On the other hand this circle lies in the union of the images of ipi ± and ip 2 ,±- 
Those functions are the unique classical solutions to 

T(i>i,±)t + F (t)y/ 1 + (ipi,±)l = 0 on R x 
ipi,±(y, 0 ) = ±\Z^o - y 2 on R x 

f T(ip 2 ,±)t +F{t)^/{ip 2 ±)l + 1 = 0 on R x 
1 ' l p 2 ,±( x } 0) = ±y/r% — x 2 on R x 

This suggests the following procedure to build AT (1) 

(Inter.) Then solve for ipi t ± and ip 2 ,± restricted to [— R,R] x 
e > 0. (2) Finally, solve for ip 3 : +. 

Some questions immediately come to mind. Criteria to decide when to move from 
(1) to the intermediate step must be chosen. Similarly, knowing which equation to 
solve within the intermediate step is a concern. The practical aspects of how a code 
reconciles the results of those steps need to be addressed carefully. We discuss all of 
these issues, and, as a result, turn the above formal idea into an efficient algorithm 
that constructs At. 

2.5. Notation. To lighten the notation, we will now work in the setting where 
n = 2. All the results discussed extend to arbitrary n. 

Continuous setting. We use the letter ip to denote functions whose image locally 
describes Ad. Suppose ip : U i-A R with ip : (y,t) ip(y,t ) = x. We introduce the 
following subsets of R 2 : 

T t := {(x,y) e R 2 : ip{y,t) = x, ( y,t ) € U} (2.12) 

See Figure 2.2 for an illustration. We distinguish between n(x, t) the two-dimensional 
outward normal to Ct at x; and n(x, t ) the three-dimensional outward normal to At 
at (x,t). 

Discrete setting. The spatial grids have fixed meshsize Ax = Ay =: h. We use 

Xi = i ■ h yj = j ■ h t k = k ■ At ( i,j , k) £ ZxZxjNU {0}} (2.13) 

to denote discrete values of space and time. We usually make no distinction between 
the continuous functions ip and their discrete approximations, except in §4. We will 
be using indices consistently, so that ipij can be understood as tp(xi,yj ) and ip k as 
ip(xi,t k )- Nevertheless, we will explicitly mention which representation is used. If a 
point p belongs to A4, then it may be described by one or more of the following three 
expressions: 


( 0,11 

{«} 


( 2 . 10 ) 


{ 0 } 


( 2 . 11 ) 


First, solve for ipz-. 
[- — e, \ + e] for some 


p) = (.fpj,yj,t k ) 


Pi — ) Pij — (xi,yj,ipij) 


(2.14) 
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Fig. 2.2. The subset To. 55 associated with ifi,— from § 2.4 appears as a black plain line. 

3. A FMM for time-dependent speeds: The t-FMM. We first address the 
problem stated in §2.1 under the following restriction: 

F = F(x, t) > 5 > 0 V (x, t) <E M 2 x [0, T] (3.1) 

Allowing the speed to depend on time yields a non-autonomous control problem. 
In [55], the author studies this min-time-from-the-boundary problem in the context 
of anisotropic front propagation. In our context, the main result of [55] may be 
formulated as follows: The value function ip for this control problem satisfies the 
following Hamilton-Jacobi-Bellman equation: 

H(Vip,ip,x) := ||VV'(x)||F(x,'i/’(x)) = 1 

The implementation of the resulting boundary-value problem: 

ll V V’(x)|| = f(x,^,(x)) < 7 °n Aq \ Cq 
ip(x) = 0 on Co 

closely mimicks that of the classical FMM. The only step that requires modifications 
is the one where a tentative value is assigned to each point in the Narrow Band. 
Following [55] this step is adjusted as follows. Let x,j = ( Xi,yj ). Without loss of 
generality, assume that Xi_iy and Xij+i are Accepted neighbours of x, ; . Consider 
a straight line lying in Quadrant II and ending at x^-, and suppose it intersects the 
line joining x,_i :? and x^j+i at the point x. See Figure 3.1. Then: x = £x ,;_ij + 
(1 — £)x ij+ i for some £ € [0,1]. Letting v = 5iij — x, we get |u| = \/? 2 + (-*- — £) 2 h- 
Associate the following value to Quadrant II: 

■0ii= min (^(x) + \/^ 2 + (1 - g) 2 (3.4) 

£€[o,i] l F(x ij ,ip(x)) J 

Proceeding similarly in the other quadrants yields the values ip i, V ; iii and ^iv- The 
tentative value assigned to ipij is then ipij = minj^i, tpu, ipm, ipiv}- Note that in two 
dimensions the minimization problem (3.4) may be solved using a direct method; see 
Appendix A. This method converges to the correct viscosity solution, and is globally 
1 st order [43, 46, 55]. Its complexity is 0(N n log N n ). In subsequent sections of this 
paper, we will refer to this modified FMM as the ‘t-FMM’. The results presented in 
this section yield Algorithm 2 given in §6. Finally, in the general case |Fj > 5 > 0, 
the PDE we wish to solve is ||V - i/ , (x)|| |F (x, ip(x)) | = 1. 


(3.2) 

(3.3) 
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Exact 

characteristic 

Approximated 

characteristic 


Fig. 3.1. If the characteristic comes from Quadrant II. 


4. A local description of the evolving front: The sideways represen¬ 
tation. An option to study the evolution of propagating curves or surfaces is to 
represent the front as a function that depends on time, e.g., y = Y(x,t) [42]. Al¬ 
though successful at describing the evolution locally, this approach fails to capture 
the global properties of the front. Nevertheless we believe that this approach can be 
used near regions where F vanishes. 

4.1. Heuristics. We first present an argument in the smooth setting. Consider 
the solution <p to IVP (2.2). Suppose <p £ C 1 (x 0 ,to) and (p(xo,to) = 0. Assume 
furthermore that <p x (ic 0 , t 0 ) ^ 0, so that the mapping is locally invertible. From 
the Implicit Function Theorem there exist open neighbourhoods (x 0 ,to) 6 V and 
U C R x [0, T], as well as a function 

i/j:U — » R , ip : (y,t) ha x = ip(y, t) , (ip(y, f), y, t) £ V (4.1) 

satisfying <p{ip{y,t),y,t) = 0 V (y,t) £ U. Taking full derivatives of cp with respect to 
y and t, and using the fact that in V, <p satisfies the LSE pointwise gives: 

{-(pxipt) + Fy/4>l + (-^xV’y) 2 = 0 <*=> -ipt ± F^jl + ipl = 0 (4.2) 

where <p x and F are evaluated at (x,y,t) = (ip(y,t),y,t). The sign used in the last 
equation depends on <p x = We let a := — sign(<^ x (xo,to))- 

Now, let ip satisfy the following Initial Value Problem: 

f ip t + aF(ip, y,t) yjl + ip\ = 0 on«n(Rx(f 0 ,T)) ^ 

\ 4>{y, t 0 ) = ipo(y) on Unfix {to}) 

where ipo is chosen such that (p(ipo(y),y,to) = 0. Then for all t £ (to,T) the set r* 
locally describes the curve at time t, i.e., r t = Ct D V. We now investigate the case 
where A4 is merely C° . For simplicity, we work with to = 0. 

Remark. Applying the same argument assuming <(> t (xo,io) 0 allows one to 

formally relate the LSE to the Eikonal equation [35]: 

(pt + F(x, y, ip)^/ (-<p t ip x ) 2 + {~My) 2 = 0 HW’M = ( 4 - 4 ) 

But since by the LSE we have a := —sign (cp t ) = sign (F(x,y,ip)), this simplifies to 

Ill’ll = \ F ( x , y , t / j )\' 
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4.2. Theory. Equation (4.3) is a Cauchy problem of the form 


f ipt + H{y 1 1, ip, ip y ) = 0 onWn(Kx (0,T)) 
l = ip 0 (y) on W n (1 x {0}) 


(4.5) 


where the Hamiltonian H : K 1 x (0, T) x K. x R 1 —> R. is defined as H(y,t,i/j,iJj y ) = 
af(i), y , t) 1 + i/jy. The function ipo is defined such that for all y £ U fl (R x {0}) we 

have {4>o(y),y) £ Co- We resort to the rich theory of viscosity solutions of Hamilton- 
Jacobi equations to study various properties of this problem [7, 8 , 14, 15, 17, 24, 28, 
49, 50]. We first address the well-posedness of the PDE. It is a simple matter to 
verify that the assumptions on H required to apply Theorem 1.1 in [48] hold in our 
context . 1 This yields 

Theorem 4.1 (Existence & Uniqueness). There exists a unique viscosity solution 
if to problem (4-5). 

We next verify that (4.5) does have the geometric interpretation advertised in the 
previous section. 

Theorem 4.2 (r t locally describes Ct). The set T t enjoys the following property: 

r t = c t nv. 

Proof. Consider IVP (2.2) again: 

/ 4> t + F\V(j)\ = 0 on K 2 x (0, T) 

\ ^(x, 0 ) = </>o(x) on R 2 x { 0 } 1 ' ' 


Since it is known that x £ Ct D V if and only if <)>(x, t) = 0, we may prove the theorem 
by showing that: x £ T t if and only if </>(x, t) = 0 . 

| =4> | We argue by contradiction. Suppose the set T = {T > t > 0 : 3x £ 
r t s.t. </>(x, t) ^ 0 } is not empty and define t* = inf T. Since (j) is continuous, T is open 
and t* fL T. Therefore, for all x* £ r t *, ^(x*,U) = 0, but for any e > 0 sufficiently 
small, there exists x e £ Tt+e such that ()>(x e ,t + e) 0. If M is differentiable 
at (x*,t*), this contradicts the argument presented in §4.1: The Implicit Function 
Theorem guarantees that the set V is open. If M is not differentiable at (x*, t*), then 
fix e and for S > 0 consider x° £ r t+e such that ||x e — x°|| <5 and M is differentiable 
at (x°,t + e). For any 6 , such a point can be found since for any T > t + e > 0 the 
singularities of r t ^_ e are subsets of measure 0 . 2 Again, the Implicit Function Theorem 
guarantees that there is a neighbourhood V of (x°,t + e) where <(>(x,t) = 0 for any 
x £ r t nV. Considering the sequence <5 rl = {^:n£N} and the corresponding 
sequence {x”})^ =1 , we arrive at the conclusion that </>(x e ,t + e) 0 contradicts the 
continuity of f>. 

| <y= | Assume that there exists (x,t) £ V such that (/>(x, t) = 0, but there is no y 
such that x = {ip(y),y) £ T*. We re-use the arguments given in the proof of | | : If 

A4 is differentiable at (x, t) then this contradicts the argument in §4.1. If A4 is not 
differentiable at (x, <), then we can find a sequence x™ £ T t converging to x such that 
</>(x", t) = 0 , and obtain the contradiction that if is not continuous. □ 

4.3. Generalizations. More generally, the above arguments can be applied to 
yield that there exists neighbourhoods (x 0 ,fo) € V and U C R x [0,T], as well as a 


1 with the exception of (H3) in [48]. However, it may be modified to get 7 R p £ R if p £ Bjy(0, P) 
for some P > 0. 

2 This follows directly from the fact that Problem (4.5) is a first order Hamilton-Jacobi equation. 
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unique function if : U —> R, if : (z,t) w = if(z,t) with (wcos(9) + z,ws'm(9) + 
z,t) £ V satisfying 

f ipt + aF (w cos (9) + z,w sin(0) + z,t) yjifl) + 1 = 0 on U D (M x (0, T)) , . 

\ if(z, 0 ) = ifo(z) on W fl (1 x { 0 }) 

in the viscosity sense. Here 9 is the polar angle of x 0 , a = —sign(xo-n(x 0 ,to)) an d V'o is 
chosen such that for all z 6 Wfl(t x { 0 }), we have (i/>o(z) cos( 0 )+z, t/>o(z) sin( 0 )+z) £ 
Cq- When 9 = 0 we recover Problem (4.3), whereas when 9 = 7 t/2 , we get that 
ip : (x, t) i— >■ y = ip(x,t) with (x,ip(x,t),t) £ V satisfies: 

f ipt + aF(x, ip, t)y/ip% + 1 = 0 on tfn(Rx (0, T)) , , 

\ ip(x, 0) = ipo(x) on U fl (R x {0}) 

in the viscosity sense. In subsequent sections, we will refer to Problems (4.3) and (4.8) 
as the yt- and xt-representations of A4, whereas Problem (4.7) will be the skewed rep¬ 
resentation. Those problems provide sideways representations of the evolving front. 
For clarity, remarks pertaining to those will usually be made for the special case of 
Problem (4.3). 

4.4. Discretization. Finite-differences schemes for problems such as (4.5) have 
been discussed [16, 18, 48]. Based on these works, we propose the following dis¬ 
cretization for Equation (4.3). In this subsection only, we will distinguish between 
the continuous function ip, and its discrete approximation which we denote as X- 
The spatial derivative \y must be computed in an upwind fashion. To this end, we 
introduce the one-sided operators 

D+ X r ■■= Xl+1 h Xl D rx r ■■= Xl ( 4 - 9 ) 

and suggest: 

X \ +1 = Xi - a- At- F(x(,yi,t r ) ■ y/l + upw(x r , l, r, a) (4.10) 

where 

upw(x r , l, n, a) := max{a, 0 } (min x r , 0 }“ + max {D ; _ x r ; 0 }“^ 

— min{ 0 , a } (max {D z + xf, 0 }“ + min {.D ; - x r , 0}“^ (4-11) 

The constant a acts as a switch and is defined as a = sign (aF(x$, yi,t r )). 

Proposition 4.3. (Convergence.) Let M be defined as the local bound on 
F, i.e., M[ = sup (a . |J/if) 6 B{p r > 2 ft ){|i r (a:,y,t)|}, where p\ = (Xi,yi,t r )- Assume that 
max {\D(~x r \i \DfiX r \} C P for all l £ L and 0 < r < R. Suppose At satisfies 

M[-At<± (4.12) 

Then the above scheme is such that x ~> if as h and At -A 0, with rate 

\\x-i>\\oo<cVAt (4.13) 

for all l, where the constant c depends on ||V’o||? ||-D^o||> the numerical Hamiltonian 
g, and RAt where 0 < r < R. 
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Proof. We proceed by showing that the scheme is monotone and consistent in 
the sense of [48]. The results then follow from Theorem 3.1 of that same paper. The 
scheme can be rewritten as 


x[ +1 = xF - A t-g ( y u t r , xh D+ X r , D( x r ) 


(4.14) 


where the numerical Hamiltonian g is easily verified to be consistent, i.e., 
9 {y,t,s, 6 , 6 ) = H(y,t,s, 6 ) V(y,t) s G R, |d| < P 
We verify monotonicity by showing that the function 


(4.15) 


G(xF-i,xF>xF+i) = XF - a- At- F(u,y h t r ) ■ y/l + upw(x r , l, r, a) (4.16) 

is a non-decreasing function of each of its argument, for fixed it, yi and t r . We only 
treat the case a > 0, since the other case is symmetric. Writing F = F(u,yi,t r ) for 
short gives 


c — aAf FyJ l +(^) 2 

if d — c < 0 , c — b < 

c — aAf F\Jl + (^jT ) 2 

if d — c > 0 , c — b > 

c — a At Fy /: L +(^) 2 + (^) 2 

if d — c < 0 , c — b> 

c — aAf F 

if d — c > 0 , c — b < 


(4.17) 


For the first case: G&, Gd > 0 are trivial to check while G c > 0 only if 


1 > \ F \ 


(At 


V h 


-1 - 


d — c 


vi+ P 2 


>M[ 


.At 


(4.18) 


Case 2 yields the same condition, whereas Case 3 gives the more restrictive one present 
in the assumption of the claim. Case 4 is trivial. □ 

Proposition 4.4. (Stability.) The above scheme is stable , provided that 

. . ( h P- 2 2 1 

\ 2 PM[ KPV 1 + 2 P 2 PS) y 1 

for some S > 0. The constant P is such that max (\D( X r \, \D( X r |} < P for all l £ L 
and 0 < r < R. 

Proof. Applying Theorem 7 of [34] to our scheme, it is possible to show that for 
h small enough, the explicit Euler map defined as 

S l At {x) = X i~ aAf • F( X i,yi,t)^/l + upw(x r , l,r, a) (4.20) 

is a strict contraction in t^. Bounding S l At ( X ) — S' At (r) from below (resp. above) 
yields the 2 nd (resp. 3 rd ) bound in (4.19). □ 

When defining ‘upw’, we implicitly assumed that both xT+i an d X'F-i were known. 
In the instance where one of those values is not known, we set X ] +i to + 00 . Indeed, no 
value can be assigned to x [ +1 since it is not possible to infer where the characteristic 
going through the point p] = { X \,yut r ) comes from. 

Remark. Assuming P = 0(l/h), we may revisit the bounds on At given in (4.19). 
The first bound is not very restrictive, even though it scales like 0(h 2 ). Indeed F « 0 
implies that M( should always be small. The bounds imposed by stability are O(h), 
which agrees with the usual CFL number of an advection problem. 
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5. Weaving the representations. Both approaches just discussed in §3 and 
§4 provide methods that locally build the manifold A4. We now address the question 
of when to use a specific representation. 

5.1. The Sign Test. Since the approach presented in §3 relies on the assumption 
that the speed is bounded away from 0, the sign of F is monitored throughout the 
algorithm. In particular, whenever a point in ( Xi,yj) is assigned a value ipij using 
the (f)-FMM, the Sign Test is performed as follows. Suppose the point Pi-ij = 
(xi-i,yj,ipi-i t j) was used in the computation of ipij. Considering the line in xyt- 
space joining the point Pi-i j and p l3 . we check the number of times d that the speed 
changes sign along this line. If d = 0, the algorithm can keep running the (t-)FMM: 
The pair (pi-i,j,Pij) is said to pass the Sign Test. If d = 1, we should change 
representation: The pair (pi-ip ,pi 3 ) fails the Sign Test. If d > 1, the grid has to be 
refined. 

5.2. Conversion of data: Interpolation. Suppose that the pair (pi~ij,Pij) 
just failed the Sign Test discussed in §5.1. Then the algorithm must change represen¬ 
tation. Without loss of generality, let us suppose that the algorithm switches from the 
xy- to the i/t-representation. This means the manifold is locally sampled by points of 
the form pi m = (a :i,y m ,t), where l € L C I and m € M C J. The yf-representation 
requires points of the form p r m = ( x,y m ,t r ), where r £ R C K. See Figure 5.1. 

5.3. Computing the outward normal. Computing the outward normal h 

accurately at each point sampling A4 is a crucial component of the algorithm. In 
regions where the level-set function cp is C 1 , we have n = ■ We use the 

Implicit Function Theorem: If ip(y,t) = x satisfies <p(ip(y,t),y,t) = 0, then (p y = 
-(pxipy and (p t = —(pxipt- Since (p x ± 0, we set n = (+s ign((p x ),ip y ,ip t ) and n = n/|n|. 
We keep track of the normal associated to each point by defining the function 

Norm : R 2 x R + — S 2 Norm (p i3 ) = n(pij) (5.1) 

5.4. The Orientation Test. Whenever a point is computed, the algorithm 
determines the orientation of the outward normal at this point. As explained in §4, 
this can be done based on the sign of fi. 3 , the time component of h. We define 

Orient^ : I 2 x 1 + —>- {—1, +1} Orients (pij) = —sign ( 713 ) (5.2) 

The algorithm requires finding which points p a b in a neighbourhood of Pij have the 
same orientation as pij. This is done using the Orientation Test. A pair (py ,p a b) 
is said to pass the Orientation Test if OrientS ( ptj ) = Orients ( p a b), and to fail it 
otherwise. 

6 . Algorithms &: Discussion. We introduce some notation before giving the 
details of the algorithms. We make use of four lists. Accepted and Narrow Band are 
lists of triplets, e.g., pij = ( Xi,yj,ipij ). Pile and Far Away are lists of coordinates, 
e.g., ( Xi,y 3 ). We define the space and time projection operators as follows: if p^ = 
(Xi,yj,ipij), then 

7r s : R 2 x R + — >R 2 n s (pij) = (x i: yj) (6.1) 

7 r t : R 2 x R + —R + 7 T t {pij) = ipij (6.2) 

The following function will be used: 

Grid : R 2 —>■ R + Grid : ( Xi,y 3 ) —> ipij (6.3) 
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t t 





Fig. 5.1. Converting the data using interpolation, (a) Some data in the xy-representation. 
The point p a p appears in red. (b) We only keep those points used for interpolation, (c) Performing 
one-dimensional interpolation line by line, we obtain data in the yt-representation. Those are the 
light green squares, (d) Using the boundary data from (c), the sideways PDE can be solved, to obtain 
the dark blue squares. By design, the domain shrinks by two points every time step. 


The set of coordinates N((xi,yj)) = {( x a ,yb ) : | {i,j) — (a, 6 )| = 1} consists of 
the nearest neighbours of ( Xi,yj ). We use Table 6.1 to define two sets of triplets: 
NeighEik((:Ei, yj)) and NeighSide(p Q ^). The first one is used to compute the value 
ipij in Algorithms 5 and 2. Similarly the second set is used in Algorithm 4, where 
the relevant component of n(p a p), the normal at p af 3 , is denoted by 17 *. We are now 
ready to present the main algorithms. 

6.1. Algorithm 1, Main loop. All steps of the main loop can be checked to 
be such that if F = F(x,y) > 5 > 0, V(x, y) € R 2 , it reduces to the classical FMM. 
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S = NeighEik((xi, y 7 -)) 

<S = NeighSide(p a/ 3 ) 

Pab — (^et) 2/6? ^Pab) 
belongs to S 
if it satisfies 

• 0a, s/6) e N((xi, Vj)) 

• Pab € Accepted 

• Grid(x a , Vb) = 'Pab 
• ( Pa/3^ Pab ) passes the Orient.Test 

• (x a ,yb) G : l e L,m G M} 

• p a i G Accepted 
• sign (rji) = sign (hi) 
where n =Norm ( p a b ) 


Table 6.1 

Definitions of two sets used in Algorithms 2, 4 and 5 


Algorithm 1 Main Loop 
1: while Narrow Band^ 0 do 


2: procedure Accept a point 

3: tl>a.p 4— min{7 Ttipij) : Pij £ Narrow Band} 

4: Grid(x a ,yp) <— IpotPi Pa/3 4- {x a ,yp, V’a/l) 

5: remove p a p from Narrow Band add p a p to Accepted 

6 : if (x a ,yp) £ Far Away then 

7: remove (x a: yp) from Far Away 


8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 


if ipap < T then 

procedure Update Pile 

for all (x a ,y b ) £ N((x a: yp)) do 
v 4— ( x a ,y b ) - ( x a ,yp ) 

if sign(t?- n(p a p)) = sign {F(j) a ,p)) or 0 then 
if Grid(x a ,y b ) = tp ab < +oo then 
Pab (x a ,y b ,ipab) 

if Orient 3 (jp a b) ^Orient 3 (p a p) then 
add (x a ,y b ) to Pile 
else if ( x a ,y b ) £ Far Away then 
add (x ai y b ) to Pile 


19: procedure Update the Narrow Band 

20: for all ( Xi,yj ) G Pile do 

21 : compute ipij and hij using Algo. 5 if F = F(x) or Algo. 2 if F = F(x, t) 

22 : p^ t- (: Xi,yj,ipij ), Norm (p^ ) 4— n{pij) 

23: remove ( Xi,yj ) from Pile. 

24: for all p G NeighEik((xj, 2 /j)) do 

25: perform the Sign Test for the pair {pij , p) 

26: if at least one pair fails the Sign Test then 

27: proceed to Algo. 3, which returns (k,l,ipki), FAIL and fi 

28: Pij 4- (Xk,yi,l/)kl), Norm (pij) 4— h, i 4- k, j 4-l 

29: Orient 3 () 4 -sign (h^Pij)) 

30: if FAIL== 0 then 

31: if 3 qij G Narrow Band with n s (qij ) = 7 r s (pij) then 

32: remove qij from Narrow Band 

33: add p^ to Narrow Band 
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Algorithm 2 Solve \Vip(x,y)\ = \ F ( x ] y ^)\ 


u± <- nt(pi±ij ) if Pi±ij S NeighEik((xj, j/j)), +oo otherwise. 
v± <- 7rt(pij±i) if Pij±i 6 NeighEik((xi,j/j)), +oo otherwise. 
0 t- [ 0 , 0 , 0 , 0 ] 
for Quadrant=l... 4 do 
if Quadrant=l then 


l/j v <- V+, Vv <- U+, T v <*- 


\F(xi,y j+1 ,^ v )\ ' 




T7U 


TVTTTT ’ 


(and similarly for other quadrants) 
if (ip v = +oo) and (if> u = +oo) then 
9 i — TOO, 

else 

9 -f- minj g [o s]{^il’ v + (1 — Qipu + \/? 2 + (1 — £) 2 (£ t v + (1 — £) r u )} 
(see Appendix A for details) 

©(Quadrant)-*— 9, 
ipij •*— min(0), Q argmin(0) 


The sideways formulations are only used when F « 0. The first procedure, ‘Accept 
a point’ is identical to the acceptance procedure in the standard FMM [41], and we 
therefore omit to discuss it. For clarity, the point accepted during this step is labelled 
as pai3 = (x a ,y^,i/jap) in the rest of the discussion. 

6.1.1. Update Pile. This step is only performed if ip a ^ is below a certain pre¬ 

defined time T to ensure that Narrow Band is eventually empty. At this stage the 
algorithm needs to decide whether a nearest neighbour (x a ,yb) of p a /3 should be put 
in Pile. To this end three criteria are used: the position, status and orientation 
of that neighbour. Simply put, line 12 has the following effect: If F{p a p) > 0 
and the considered neighbour lies inside the curve then the pair ( x a ,yb ) is 

not added to Pile. Next the status of this nearest neighbour is considered. If the 
pair (x a ,yb) was traversed by the curve in the past, then it is only added to Pile if 
Pab '■= ( x a ,yb,Grid(x a ,yb )) and p a p have different orientations (lines 13-16). In¬ 
deed a point in the plane can only be traversed twice if the speed has changed sign 
in the meantime. If ( x a ,yb ) is still in Far Away, then it is automatically added to 
Pile (lines 17-18). 

Remark. The presence of the ‘if i/j a /3 < T' in line 8 is in contrast with the 
standard FMM, where it is proved that since F > 5 > 0, all characteristics exit 
the domain in finite time. In this context, the size of the computational domain 
determines T. 

6.1.2. Update the Narrow Band. This procedure assigns tentative values to 
the points in Pile using either the standard FMM (see Appendix B) or Algorithm 2, 
depending on the domain of F. Since i/j only solves the Eikonal equation in regions 
where |F| > S > 0, the first lines of those algorithms ensure that the points involved in 
the computation of ijiij all lie in one such region. The steps outlined in lines 24-28 
represent the main modification to the standard FMM algorithm. The Sign Test is 
performed to check if the value returned by Algorithm 5 or 2 is valid. If it is not, 
then Algorithm 3 is called. Using a sideways representation, it attempts to return a 
point {xk,yi,ipki) £ C-iAfcj- If h manages to do so, note that as explained in §6.2.4, the 
triplet returned may not be (xi,yj,ipij), which is why i and j are relabelled in line 
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28. As in the standard FMM, if there already is a point in Narrow Band with the 
same spatial coordinates (x;, yj ), then it is automatically removed from that list. The 
triplet (xi 7 yj,ipij) is added to Narrow Band. In the event where Algorithm 3 fails, no 
new point is added to Narrow Band. 

6.2. Algorithm 3, Sideways representation. This algorithm is called by the 
main loop when the speed F is close to 0. 

6.2.1. Determine representation. In order to work locally, the first step of 
this procedure defines a square of side length at most 2 sh for some s £ N as the new 
computational grid. Then the representation is chosen based on the normal at p a p. 

6.2.2. Initialization. This is the step where data are converted, as was men¬ 

tioned in §5.2. The set NeighSide^c^) is found; This ensures that the orientation of 
the points used next is compatible with the current representation. We take time to 
explain what we mean in line 12 in details. It is ideal to build the sideways grid in 
such a way that the triplet p a p is represented exactly on this grid, i.e., For example, 
if data are being converted to the yf-representation, then there should be l € L and 
f € R such that (ypt r ) = (yp,ifap)- The function if\ : (y,t) '-A x then satisfies 
ipi(yj,t r ) = x a , and (x a7 yp,if a p) = (-01 , t r ), yj, t r ). This avoids rediscovering the 

point p a p in the procedure ‘Get (ccfe, 2 /z, )’ discussed in §6.2.4. Assigning values to 
the sideways grid in line 13 is an interpolation problem. See Figure 5.1 (c). 

6.2.3. Main loop. The sideways PDE can now be solved. As mentioned in §4.4, 
if either or V'z+i 1 are set to +oo, then Algorithm 4 sets i/jJ to + 00 . As depicted 
on Figure 5.1 (d), this has the effect of shrinking the size of the set where the PDE 
is solved: At most s time steps can be taken before all the boundary information 
available has been used up. When the speed depends on time, we believe that using 
adaptive time stepping increases the success rate of Algorithm 3. We pick a small At 
as long as the speed has not changed sign. This makes the scheme more accurate, 
thereby increasing the chances of assigning a value to ( Xi,yj ). Once F changes sign, 
a large At is chosen to increase the likelihood of assigning a value to ( x a ,yp ). 

6.2.4. Get ( Xk,yi 7 Pki )• Deciding which value is returned by the algorithm is 
delicate and may be summarized as follows: By default, the algorithm always tries to 
assign a value to the pair in the Narrow Band (lines 22-25). If this is not possible, 
then it tries to assign a new value to the pair (x a , yp) = Tf s (p a p) (lines 26-29). If this 
cannot be done either, then this representation failed. The algorithm must attempt 
using another representation which is chosen based on the ones already attempted. 
When the 1 st attempt fails. Suppose the xt-representation failed, then the algorithm 
attempts to use the yf-representation. When the 2 nd attempt fails. Then the scheme 
resorts to the skewed representation. When the 3 rd attempt fails. If the skewed 
representation also fails, then Algorithm 3 fails entirely. Note that this is expected to 
happen if ( Xi,yj ) and ( x a ,yp) Ct for any t £ ( p a p,T ). See Example 2 in §8. 

Remark. In practice, after each iteration of the for loop line 17, we check if 
either ( Xi,yj ) or ( x a ,yp ) has been traversed by the curve. If not, then the for loop 
keeps going. 
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Algorithm 3 Sideways representation 


procedure Determine representation 

s € N is picked, v t— (— s, —s+l,...,s—l,s) 

L 4 — L n I , AT 4 — AT n J 

if \rii{p a p)\ > \n 2 (p a p)\ then 

use yt-representation: z 4— x, a 4 -sign(hi(p aj g)) 

else 

use xf-representation: z 4— y, a 4 -sign(n, 2 (p a ^)) 


L 4^- i + v, AT 4— j + v 


8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 


Attempt 4— 1 
while Attempt> 0 do 

procedure Initialization 
get NeighSide(p a/3 ) 

the sideways grid ( zi,t r ), l G L, r S R is built 

Grid(zi,t r ) t— ipl using interpolation and NeighSide(p a| g) where possible. 
Grid(zi : t r ) 4 —t-oo where interpolation cannot be used, 
procedure Main loop 
if a ^ 0 then 

for n = 1 : R max do 
At is determined 
for l = 2 ; L max 1 do 

compute 'ipj using Algo. 4. 

procedure Get (x k , yi,4>ki) 

if ( Xi,yj ) is traversed by the curve then 
ipij is computed using interpolation 
4>ki ipij, x k 4- Xi, yi 4- yj, h(ipkl) is computed 
Attempt 4— 0, FAIL 4— 0 
else if (x ai yp) is traversed by the curve then 
ip a p is computed using interpolation 
4>ki 4- Ip a p, x k 4- x a , yi 4- yp, h{4>ki) is computed 
Attempt 4— 0, FAIL 4— 0 
else This sideways representation failed, 
if Attempt =1 then 

if in xt-representation then 

use i/f-representation: z 4— x, a 4 -sign(ni(p a ^)) 

if in yf-representation then 

use xt-representation: z 4— y, a 4 -sign(n 2 (p a / 3 )) 

Attempt = Attempt +1 
else if Attempt=2 then 

use the skewed representation: z 4— w, a 4 -sign((x a , yp)-n(p a p)) 

Attempt = Attempt +1 
else if Attempt=3 then 

Point is not reached before T. ipki < —boo, x k 4 —boo, yi 4 —boo 
Attempt 4— 0, FAIL 4— 1 






18 


A.TCHENG, J.-C.NAVE 


Algorithm 4 Solve ip t + aF(iJj, y , t)^J 1 + if } 2 = 0 

if (ipili < +oo) & (ipi+-1 < +oo) then 
a <- sign (aF(il>l~ 1 ,y h t r - 1 )) 

Vi t- Vi~ X ~ a- At ■ F('t/jl~ 1 ,yi,t r ~ 1 ) • ^/l + upw (i/j r ~ 1 ,l,r, a) 

else 

< - 1-00 


6.3. General remarks. 

6.3.1. Data structure. One of the main differences with the standard FMM is 
the way we keep track of the various properties associated to each point. The fact that 
a point (x a ,yp) on the plane may be traversed by the curve more than once requires a 
slightly richer data structure. For example, the functions Norm and Orient 3 have to 
be defined over triplets rather than over R 2 . On the other hand, the lists Pile and 
Far Away still consist of coordinates. Note that when the code ends, Narrow Band is 
empty whereas Far Away may still contain points. The Accepted list may contain 
multiple triplets sharing the same spatial coordinates. In order to keep track of what 
the most ‘up-to-date’ value associated with ( x ai yb ) is, we make use of Grid. Indeed, 
this function enjoys the following property: If there are distinct points Pij, qij G 
Accepted such that 7r s (py) = ir s (qij), then Grid(xi,yj) = m&x{Trt(pij),nt(qij)}. 
Viewed as a set, Grid(7 r s (Accepted)) is the upper semi-continuous envelope of At. 

6.3.2. Recovering the curve from Ai. The set Accepted provides a discrete 
sampling of Ai. Using this point cloud, and possibly the normal h to Ai at each point, 
a continuous representation of Ai can be obtained. See for example [6, 10, 38, 27, 30], 
and [56]. Given a time t € (0, T), a contouring algorithm can then be used to find C* 
(see [30]). 

6.3.3. Resolution. By construction, the density of points sampling Ai is ex¬ 
pected to be lower in regions where F « 0. A remedy to this situation is to also 
record the points computed in the sideways representations. 

7. Complexity of the method. We derive some estimates for the computa¬ 
tional time of the method when n = 2, i.e., two spatial dimensions. Consider a spatial 
grid of N 2 points with meshsize h. Let At ~ h, and define N* to be the number 
of gridpoints traversed by C* when 0 < t < T. (i.e., if a given gridpoint ( Xi,yj) is 
traversed twice, say at times t\ and t 2 where 0 < t\ < t 2 < T, then this contributes 
+2 to N*.) By construction, the computational time depends on the size of the set 
Fm := F G Ai. Indeed, Algorithm 3 is only called when Algorithm 1 fails, which 
occurs whenever an accepted point computed by Algorithm 1 is within a spatial dis¬ 
tance h of Fm ■ Let the number of points computed by Algorithm 3 be N. Since the 
complexity of Algorithm 1 is well-known [41], let us focus on estimating the complex¬ 
ity of a single call to Algorithm 3. On the square of side 2s, the Narrow Band forms a 
one-dimensional subset. Using interpolation to convert the points in a neighbourhood 
of this set takes 0{s) operations. Algorithm 4 makes at most s 2 operations. Those 
two steps are performed at most three times. We formally argue that the parameters 
of the algorithm can be chosen such that this worse case complexity is not achieved. 
The procedure mentioned in the remark of §6.2.4 can be used to prevent Algorithm 
4 from making unnecessary computations. In §6.2.3, we explain how using adaptive 
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time-stepping increases the success rate of Algorithm 3. Moreover, as N increases, 
the time distance between the accepted point computed by Algorithm 1 and Fm de¬ 
creases, which in turn makes Algorithm 3 more successful on average. Altogether, this 
suggests that the number of attempts taken by Algorithm 3 tends to one for almost all 
points; this is confirmed by the examples presented in the next section. As a result, 
the complexity of Algorithm 3 tends to O(s) for large N. Given the assumption that 
F is analytic, we expect N* — N = 0(N 2 ) and N = O(N). In practice, the number 
of points in the local grid s can be chosen as kN for k -C 1. The overall complexity 
can therefore be estimated as: 

0(N 2 \og{N 2 )) + 0(N)xO(kN)=0(N 2 \og{N 2 )) + 0{kN 2 ) (7.1) 

(t)-FMM augmented part 

Note that in the instance where F = 0, we recover the usual complexity of the FMM, 
namely 0(N 2 \og(N 2 )). 

8 . Numerical Tests. In this section, we illustrate how the method works with a 
variety of examples. We first discuss the methodology used to assess the convergence 
of the algorithms, and briefly summarize which features and results are expected. We 
then present the examples. More details are provided in Appendix C. 

8.1. Error measurement. To assess the convergence of our algorithm, we com¬ 
pute the error associated to each point pij returned by our scheme. 

Method 1: Eij. Suppose that an exact solution to the Level-Set Equation (2.2), 
cf>(x,y,t) > 0 is known, with the property that |V<^| = 1 for all t. Then evaluating <f> 
at = (xi,yj,ipij) returns the distance to the curve C^ ir We define Ejj = \<f>(pij)\. 
This method is used for all examples except Example 4 when F < 0. 

Method 2: Gij . If an exact solution is not available, we get a numerical solution 
accurate enough to be considered exact. To this end, the Level-Set Equation is solved 
on a very fine grid using 2 nd order stencils in space, and RK2 in time. At each 
time step, the zero-contour of <f> is found and sampled. The resulting list of points 
B provides a discrete approximation of A4. The error associated to p l; j is defined 
as the smallest three-dimensional distance to this exact cloud of points, i.e., Gij = 
mm qe i 3 {\pij — g|}. This method is used for Example 4, when F < 0. 

8.2. Tests performed. 

Accuracy of Algorithm 4- In §4.4, it is mentioned that the sideways method we 
propose converges with at least 0(h 1 ^ 2 ) accuracy. To verify this, we pick a domain 
U , initialize say x = ^(z/mi^ 0 ) with exact data for some initial time t°, and run 
Algorithm 4 for different gridsizes. The result is a subset of M, encoded as a list 
of points of the form p r m = ('iffn, y m , t r ). An error is associated to each point p r m 
such that iffn < oo using either Method 1 or 2, i.e., E r m = Inexact (ym,t r ) — tffnl or 
G r m = minqgBdp^ — g)}. A two-dimensional L\ norm is then used to report the 
results in Figure 8.1, e.g., L x = h 2 ■ Erefl Bin- 

Accuracy of the full scheme. When testing the accuracy of the full scheme, we 
distinguish between different regions of the resulting set Accepted. When study¬ 
ing a region computed by the (t)-FMM, a two-dimensional Li norm is used: Li = 
h 2 • E ie / EjejEij- Note that our assumptions on F imply that the points com¬ 
puted using the sideways representations form one-dimensional sets of R 2 x [ 0 ,T]. 
Consequently, a one-dimensional L\ norm is used to study those points: L\ = h ■ 
J2i£l J2jejEij. The global error (computed using all the points in Accepted) is a 
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two-dimensional Li norm. It may be interpreted as an approximation of the volume 
enclosed by the exact and the approximated surfaces. 

We report the L a0 error qualitatively, through the black & white representations 
of the set Accepted. Those figures are obtained by computing the relative error at 
each point, i.e., if = max^gj, then = Eij/L ^; and then shading the 

point accordingly: The darker a point, the larger its relative error . 

8.3. Expectations. By assumption, as h —> 0, the 1st order (f)-FMM scheme is 
used almost everywhere. This should reflect in the global error: It should follow the 
same trend as the (t)-FMM. Moreover, we expect the call to Algorithm 3 to increase 
the constant of convergence. A question that we address is the extent to which this 
degrades the local and global accuracy. We investigate the behaviour of the scheme 
in the presence of shocks & rarefactions in Example 4, as well as in §9. 

In all examples but the fourth one, the initial curve Co is the circle centred at the 
origin, with radius ro = 1/4. In all tests, data are initialized with exact values. 

8.4. Example 1 : F = F(t) = 1 — e 10t_1 . The main purpose of this example is 
to illustrate the basic ideas at play in the method. The speed is such that we expect 
the circle to first expand up to time t = 0.1 and then contract until it collapses to 
the origin. We first assess the order of convergence of the method for the sideways 
representation. The results reported on Figure 8.1 clearly indicate that it is 0(h). 
This is higher than the (^(/i 1 / 2 ) rate that was predicted in §4.4. When the entire 
code is run, the set of Accepted points is presented on Figure 8.2 (a)-(b). One¬ 
dimensional optimization is used for those points traversed by a characteristic that 
is almost aligned with one of the spatial axis. We note that the sideways points 
are computed in the yt- (resp. :rt-)representation when n aligns better with the x- 
(resp. y-)axis. As expected, the sampling of the surface is sparser near the plane 
t = 0.1. Remark that in this example, none of the sideways points were computed in 
the skewed representation. The global convergence results are presented in Figure 8.2, 
(d). We distinguish between the bottom part of the surface, the top part, and those 
points computed using the sideways representation. On the one hand, the results 
pertaining to the bottom part allow us to conclude that the t-FMM is O(h), as 
predicted in §3. On the other hand, we can study the effect of the call to Algorithm 3 
on the behaviour of the scheme. Indeed, although the f-FMM also converges with O(h) 
when used to build the top part of the surface, it does so with a larger constant. We 
conclude that changing representation does deteriorate the accuracy of the sampling 
but only to a mild extent. To gain a better understanding of where the loss of 
accuracy from the bottom to the top part stems from, the relative L^ error e+j 
associated to each point can be viewed on Figure 8.2 (c). Those points computed 
using one-dimensional optimization in the f-FMM, just after t = 0.1 bear the largest 
errors. Two reasons can explain this: Some of those points are clearly traversed by 
characteristics that are not aligned with the spatial axes. Nevertheless, the scheme 
resorts to one-dimensional optimization to assign them values, for lack of a better 
method. Indeed, when those points are put in Pile, there are not enough neighbours 
with negative orientation available to use two-dimensional optimization. We also 
suspect the constant of convergence of the t-FMM to depend on S where |Fj > <5 > 0. 
In practice, this method is found to perform poorly when using points pij such that 

F{pij) ~ 0. 

Remark. Those outliers do not degrade the accuracy of the method, even locally. 
This is because by design, Fast Marching Methods assign values to the points in 
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Fig. 8.1. Convergence results for the sideways scheme, using (left) Method 1 (right) Method 2. 


Pile using only those neighbours with a smaller value. As a result, those outliers are 
not used in any of the calculations of the values of their neighbours. In practice, it is 
found that they eventually become isolated points of the Narrow Band before getting 
accepted. 

8.5. Example 2: F = F(x) = x. The given speed is such that the curve remains 
a circle whose radius grows while its center shifts to the right. Our method adequately 
handles this case as a single problem, although the speed changes sign across the y- 
axis. As expected, Algorithm 3 fails near the points (0, 0.25) and (0, —0.25), as shown 
on Figure 8.3 (a). The sideways scheme was tested both in the xt- and the yt-charts, 
and was found to be 1 st order in each case (Figure 8.1). The results for the full scheme 
show that it converges with O(h) accuracy everywhere (Figure 8.3 (b)). Let us bring 
up that a bi-directional FMM was proposed in [13] to solve a related problem. 

8.6. Example 3: F = F(x,y,t). (See Appendix C.5 for details about F.) This 
example differs significantly from the previous ones in that the set F no longer consists 
of planes. The exact solution C t is a circle that only grows at first, and then starts 
moving in the positive ^-direction. Our method is observed to perform very well on 
this example; We present the resulting surface and the first order convergence results 
on Figures 8.1 & 8.4. 

8.7. Example 4: Two merging circles. This example tests the ability of the 
scheme to capture topological changes. The initial codimension-one manifold consists 
of two disjoint circles of radius ro = 1/4, with centres at (—.3,0) and (.3,0). The 
speed is such that the circles first expand, until they touch and merge. Then the 
speed changes sign, which makes the curve shrink until it pinches off and splits into 
two distinct curves. The set Accepted is presented in Figure 8.5 (a).The accuracy 
of the sideways scheme is investigated on a domain that comprises the shock when 
F > 0, and the rarefaction when F < 0. First order convergence is obtained in each 
case (Figure 8.1). The full scheme also shows 1 st order convergence (Figure 8.5 (b)). 
The convergence of the sideways points and the top part is a little shy of first order, 
but this can be attributed to the measurement method. Those results demonstrate 
how robust the overall scheme is. Note that a similar example was tackled in [9], with 
a speed F that depended linearly on time. 
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Legend for (a),(b): 



Initialization 


★ 


i-FMM, ID 


t-FMM, 2D 



yt -representation 
xt -representation 




Fig. 8.2. Example 1: (a) - (b) Different perpectives of the set Accepted, (c) Accepted is 
featured. The relative error etj determines the shade of each point, (d) Convergence results. 



(a) 



Log(h) 

(b) 


Fig. 8.3. Example 2: (a) The set Accepted (b) Convergence results. 
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Fig. 8.4. Example 3: (a) The set Accepted (b) Convergence results. 




(a) 


(b) 


Legend for (a): t -fi 

Initialization 


•FMM, ID 
£-FMM, 2D 


□ 


yt -representation 
x t -representation 


Fig. 8.5. Example 4 : ( a ) The set Accepted, (b) Convergence results. 


9. Discussion. In the light of the examples presented in the previous section, 
we address the limitations, weaknesses and advantages of the algorithm. 

We illustrate one of the main limitation of the scheme with an ultimate example. 
The speed is chosen such that the initial circle immediately develops a kink along 
the a;-axis at time t = 0. Its subsequent shape resembles that of an almond slowly 
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y 


x 


x 


Fig. 9.1. A4 for the almond example. 


The shock appears as a red plain line. 




(a) ' (b) 


Fig. 9.2. The almond example: The set Accepted (a) side view, (b) viewed from above. 




turning in the counterclockwise direction while expanding. The sign of the speed 
changes, forcing the curve to contract while retaining its slanted shape. See Appendix 
C.7 for details. The most prominent feature of this example is that, as is depicted 
on Figure 9.1, the shock is not a straight line. Remark that the speed F does not 
satisfy the assumptions of this paper outlined in §2.2: It is only a C° function of 
R 2 x [0,T]. The surface that results from running the algorithm at high resolution is 
shown on Figure 9.2. The shock is clearly visible, and has the expected figure-eight 
shape. Nevertheless some points ‘escape’ through the shock when the speed changes 
sign, and start out two new fronts that keep on expanding. The problem stems from 
the procedure ‘Update Pile’ in Algorithm 1. In order to decide which points go in 
Pile, the code distinguishes between the inside and the outside of the curve using 
the normal n (cf. line 12 of Algorithm 1). Consider what happens along the shock, 
where n has a discontinuity. So long as the expansion is outwards, this does not 
cause problems. But when the direction of propagation reverses, some points that 
should stay in Far Away are moved into Pile. A possible remedy to this issue is to 
approximate the normal cone along the shock. This additional information could be 
used as an updating criterion. 
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y 


• in NeighSide(p aj g) 
o Accepted points 


Fig. 9.3. Illustration of what happens if s is chosen too large. Data need to be converted to the 
xt-representation, but the set NeighSide(p a p) of neighbours of the black point ( x a ,ya ) consists of 
two connected components. 



On a much more general note, the gluing mechanism between the two formalisms 
heavily relies on an accurate computation of the normal. In practice, we found that 
the algorithm is rather sensitive to the accuracy of this quantity. Another weakness 
of the method is that, as it stands, Algorithm 3 may fail when it is not supposed to. 
i.e., Even though ( Xi,yj) or (x a ,yp) belongs to Ct for some t £ (0 ,T), the algorithm 
does not assign any value to either of those coordinates. Two situations make such a 
scenario possible: (1) the time steps taken are too small, or (2) too little information 
obtained from interpolation is available. Recall from Proposition 4.3 that the CFL 
condition prevents large At. Case (2) can occur if s £ N, the number of points 
in the local grid in Algorithm 3 is too small. However, if s is large, Algorithm 3 
may not be able to carry out the step outlined in line 13. This happens if the 
points in NeighSide^c^) sample more than one connected component of the set {p £ 
A4 : 7 r s (p) £ [xi- s ,Xi+ s ] x [yj- s , 7/j+ s ]}. See Figure 9.3 for an illustration. However, 
choosing s systematically so as to prevent this situation seems difficult. Ultimately 
h and s depend on measurable quantities such as the Lipschitz constants of F and 
its derivatives, as well as the local curvatures of T t . Nevertheless the way those 
parameters are intertwined and should be chosen is a question we wish to address in 
future work. 

The fact that our method is a rather mild modification of the standard FMM 
has obvious benefits. As featured in all the examples, the sideways representations 
need only be used to compute a relatively small number of points sampling M.. This 
allows us to safely predict that the computational complexity of the algorithm is lower 
than that of pre-existing algorithms used to tackle this problem, such as the LSM or 
the GFMM. Nonetheless, it is hard at this point to make more precise complexity 
statements. 

10. Conclusion. Our aim was to devise an algorithm with low complexity able 
to describe the non-linear evolution of codimension one manifolds subject to a space- 
& time-dependent speed function that changes sign. To this end, we illustrated how 
pre-existing methods can be combined to achieve this goal. The fact that we always 
dealt with explicit representations of the manifold implied that the dimensionality of 
the problem was never raised. The resulting algorithm was found to have a global 
truncation error of 0(h). We tested it against a number of examples, some of which 
cannot be found in the current literature. 
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The algorithm is found to be robust and accurate in all the tests presented. 
Regarding the complexity of the method, a legitimate concern is to clearly quantify 
how the success rate of Algorithm 3 depends on the various parameters involved, as 
well as the speed function F and the manifold AT Once this is done, more precise 
statements about the runtime of the algorithm can be made and tested. 

Overall, the present work thoroughly introduces a new algorithm, along with 
proofs of convergence and stability, as well as sturdy numerical results. We believe 
that the main idea on which it relies - i.e., to change representation based on the speed 
function F - may be extended and improved in many ways that shall be explored. 

Appendix A. A direct method to compute i/ii in the t-FMM, in 2D. 

We provide a direct method for solving the minimization problem appearing in 
Equation (3.4), in two dimensions. Introducing r{y) = h ^( y ))\ > we first use linear 

interpolation to simplify the quantity we wish to minimize: 

^) + VeT (TV^_A_ = «,(*) +v? + (i-0Mx, 

« + (1 - OVKxij+i) + + (€r(xi-u) + (1 - £Mxi,j+i)) 

=: m (A.l) 

Minimizing / over £ £ (0,1) amounts to finding the roots of 0 = C4A 4 + C3A 3 + C2A 2 + 
ci A + co where A £ (0,1) is such that /'(A) = 0. This quartic can be solved either 
directly with closed formulas, or with Newton’s method — we use the latter. For 
each root r 7 ; £ (0,1) the corresponding value of '0 is computed as ipn.r, = /(rj). If 
iftu >ri < ipfa-ij) or i/’n.n < then i/’n.n is discarded. Values arising from 

minimization in one dimension are also computed as ^11,0 = i/Kxij+i) +r(xij+i) and 
■011,1 = + t(xi_ ij). The global minimum is found by comparing all those 

values. 

Appendix B. Algorithm 5, standard FMM. We revisit the standard Fast 
Marching Method algorithm, using some of the notation we have introduced. 

Algorithm 5 Solve \Wip{x,y)\ = p^yy 

u± £- TTt{pi±ij) if Pi±ij £ NeighEik((a:i, yj)), +00 otherwise. 
v± •£- n t (pij± i) if Pij±i £ NeighEik((a;i, j/j)), +00 otherwise. 
u <— min(u_,u+), v t— min(r_, n+), 
if max(n, v) — min(it, v) < ^ F ^ x h y .^ then 

tv = I + «) + \J 2 {nk^)) 

else 

ipij =mm(u,v)+ {F(x h . tyj)i 


Appendix C. Implementation details for the examples. 

C.l. Solvers used. We give some details about the examples presented in §8. 
All tests were performed using Matlab® [1], In particular, finding the minimum 
value in the Narrow Band is done using the command min . 

Whenever a value ipij is computed by the (f)-FMM, the normal fiij is approxi¬ 
mated using the one-sided derivatives involving the points used in the computation 
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of tpij. For example: if two-dimensional optimization was used in Quadrant III to 
obtain t/v,-, then 


ij 'Ipi—l.j tdj 'Ipi.j — l 

h ’ h~ 


Orient 3 (jpij) ) and n{pij) = — (C.l) 


Within Algorithm 3, we approximate the normal as follows. For clarity, say the points 
Pj = (V’j an d p k ~ l = (r/^A 1 , yj, t fe_1 ) computed in the yt-representation with 

ar-orientation a were used to obtain pij = ( Xi,yj,tpij). Then 


v = —a, a 




2 h 


-,a 


ip 1 * - ip 


k -1 ' 


dt 


and n(pij) = (C.2) 


Note that this is not an approximation of the true normal at Pij , which is (—a, —cj) x ipy, 
—4> x il>t)- However, the only two salient information we need from h are: the sign of 

fi 3 and the direction of n. The two-dimensional normal is simply obtained from h as 

- _ (m,»2) 

|(rti,ra2)| ‘ 

C.2. Choice of parameters. In all examples, the number of points in each 
dimension is N + 1, and the spatial grid spacings are even: h = dx = dy. The size of 
the local grid in Algorithm 3 was set to be s = [yj- As discussed in §6.2.3, we use 
adaptive time-stepping, in those examples where F depends on time. In the fine part, 
before the time where F = 0, we set At = rq/i. Passed that time, we let At = r 2 h. 
To assess the convergence of the sideways methods, a yt -grid with spacings h and 
At = h/2 was built.Remark that the exact normal h was assigned to the points as 
they were accepted in all the examples, except Example 1 where it was computed as 
explained in §C.l. 

C.3. Example 1. The exact solution to the Level-Set Equation is <j>(x,y,t) = 

y/x 2 + y 2 — R(t) where R{t) = (r 0 — e "q " 1 + t'j . Domain: [—.321, .319] 2 . Tp = 0.3. 

xt- and yt- rep.: rq = 1/3, r 2 = 2. Skewed rep.: rq = r 2 = 1. Domain for convergence 
of Algo. 4: (y,t) e [-0.25,0.25] x [0,0.3], 

C.4. Example 2. The signed distance function to the curve C t is given as 

(j)(x,y,t) = y/(x — x c (t)) 2 + y 2 — r(t) where x c (t) = r 0 sinht and r(t) = r 0 coshf. 

Note that (j) does not solve the Level-Set Equation. Domain: [—1.01, 0.99] 2 . Tp = 1. 
xt- and yt- rep.: rq = 1/3, r 2 = 2. Skewed rep.: rq = 1/3, r 2 = 5. Domain for 
convergence of Algo. 4: (y,t) € [—0.25,0.25] x [0,1] and (x,t) £ [—0.25,0.25] x [0,1]. 

C.5. Example 3. The exact solution to the Level-Set Equation is <j>(x,y,t) = 
y/(x — gt ) 2 + y 2 — (ro + ct) where b = 10, c = 1/2 and g[t) = arctan (b(t — 0.5)) + §. 
The speed is 


(x - gt)(g't + g) ^ 

V(x-gt ) 2 + y 2 



c 

(x-nt)n 

y/{x-Trt) 2 +y 2 


for t small 
for t large 


(C.3) 


We expect the circle to first expand (when t is small), and then expand while moving 
to the right with speed n (when t is large). Domain: [—1.51,+1.49] 2 . Tp = 0.5. xt- 
and yt- rep.: rq = 1/3, r 2 = 2. Skewed rep.: rq = 1/3, r 2 = 5. Domain for convergence 
of Algo. 4: (y,t) G [-0.25,0.25] x [0,0.5], 
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C. 6 . Example 4. The set C 0 consists of two disjoint circles of radius ro = 0.25, 
with centres at (—0.3,0) and (0.3,0). The speed is F = 1 — e 2t_1 . The circles touch 
along the y -axis when t « 0.08. When t < 0.5 the exact solution to the Level- 
Set Equation is cj>(x,y,t) = min ( yj(x + 0.3 ) 2 + y 2 — R(t), y/(x — 0.3 ) 2 + y 2 — i?(i)j 

where R(t) = r 0 — e2 *~ 1 + t. Domain: [—1.5 + 0. Ole, +1.5 + 0. Ole] 2 . T F = 1.2. xt- and 
yt- rep.: 7 q = 1/3, r 2 = 2. Skewed rep.: rq = 1/3, r 2 = 5. Domain for convergence of 
Algo. 4: (y,t) <E [-0.5,0.5] x [0.2, 0.5] and (y,t) e [-0.5, 0.5] x [0.5, .52], 

C.7. The Almond example. The exact solution to the Level-Set Equation is 

2 /, t) = (y/x 2 + y 2 -r 0 + - -- - t(l + C)\ + %== (C.4) 

V ce / Vl + t 2 

=-4>(x,y,t)+g(x,y,t ) (C.5) 

The constants are set to be: ro = 1/4, c = 1, and C = .65. The function <f> is made 
up of two parts: <j> is qualitatively the same as in Example 1. Domain: [—0.5, 0.5] 2 . 
Tp = 1.9. xt- and yt- rep.: n = 1/3, r 2 = 2. Skewed rep.: rq = 1/2, r 2 = 6 . 
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